


typename<-"stationary"

setwd("./output")

obsdraws<- as.data.table(read.dta13(paste0("psid_obs",typename,".dta")))
obsdraws_cons<-as.data.table(read.dta13(paste0("psid_obs_cons",typename,".dta")))
#NEED TO CORRECT DATAMOMENTS AND INVVARCOV:
temp<-read.table(paste0("wagemomentsstationary_jump.txt"),row.names=1,fill=TRUE)
temp<-temp[rownames(temp)%in%c("prodtype_1","prodtype_2","prodtype_3","age","age2","mod","sev","old","mod_old","sev_old"),]
wageregmoments<-temp[,1]
wageregse<-as.numeric(temp[,2])
temp<-read.table(paste0("spwagemomentsstationary_jump.txt"),row.names=1,fill=TRUE)
temp<-temp[rownames(temp)%in%c("prodtype_1","prodtype_2","prodtype_3","sp_age","sp_age2","sp_mod","sp_old","sp_mod_old"),]
spwageregmoments<-temp[,1]
spwageregse<-as.numeric(temp[,2])

temp<-read.table(paste0("wagevarsstationary_jump.txt"),row.names=1,fill=TRUE)
wagevarmoments<-temp[,1]
wagevarse<-temp[,2]


temp<-as.data.table(read.table("consmoments_all.txt"))
temp$V2<-as.numeric(as.character(temp$V2))
temp$V3<-as.numeric(as.character(temp$V3))
#dropping PT work for spouses:
temp<-temp[!grepl("sp_part_part",temp$V1),]
temp<-temp[!grepl("part_part_",temp$V1)]
consmoments_all<-temp[,V2] 
consse_all<-temp[,V3] 
temp<-as.data.table(read.table("consmoments_single.txt"))
temp<-temp[!grepl("part_part_",temp$V1)]
temp$V2<-as.numeric(as.character(temp$V2))
temp$V3<-as.numeric(as.character(temp$V3))
consmoments_single<-temp[,V2]
consse_single<-temp[,V3]  
temp<-as.data.table(read.table("consmoments_married.txt"))
#dropping PT work for spouses:
temp<-temp[!grepl("sp_part_part",temp$V1),]
temp$V2<-as.numeric(as.character(temp$V2))
temp$V3<-as.numeric(as.character(temp$V3))
consmoments_married<-temp[,V2] 
consse_married<-temp[,V3] 


temp<-as.data.table(read.table("consmoments_simple_all.txt"))
temp$V2<-as.numeric(as.character(temp$V2))
temp$V3<-as.numeric(as.character(temp$V3))
#dropping PT work for spouses:
temp<-temp[!grepl("sp_part_part",temp$V1),]
temp<-temp[!grepl("part_part",temp$V1)]
consmoments_simple_all<-temp[,V2] 
consse_simple_all<-temp[,V3] 


##dropping PTwork of head and partners:
#consmoments_all<-consmoments_all[-c(4,11,14,15,18,20)]
#consse_all<-consse_all[-c(4,11,14,15,18,20)]
temp<-read.table("workmoments_agg.txt")
temp<-temp[!grepl("E=1/",rownames(temp)),]
workmoments_simple<-temp[,1]
workse_simple<-temp[,2] 


temp<-read.table("workmoments_1.txt")
temp<-temp[!grepl("E=1/",rownames(temp)),]
workmoments_mar<-temp[,1]
workse_mar<-temp[,2] 
#Dropping PT work:
#workmoments_mar<-temp[1:6,1]
#workse_mar<-temp[1:6,2] 
temp<-read.table("workmoments_0.txt")
temp<-temp[!grepl("E=1/",rownames(temp)),]
workmoments_sing<-temp[,1]
workse_sing<-temp[,2] 
# workmoments_sing<-temp[1:6,1]
# workse_sing<-temp[1:6,2] 
temp<-read.table("spworkmoments_1.txt")
#dropping PT work for spouses:
temp<-temp[!grepl("E=1/",rownames(temp)),]
spworkmoments<-temp[,1]
spworkse<-temp[,2]
#spworkmoments<-temp[1:4,1]
#spworkse<-temp[1:4,2]
temp<-read.table("stockDImoments_agg.txt")
stockDImoments_simple<-temp$EST[]
DIse_simple <- temp$SE[]

temp<-read.table("stockDImoments_work.txt")
temp<-read.table("stockDImoments.txt")
stockDImoments<-temp$EST[]
DIse <- temp$SE[]
#stockDImoments<-temp$EST[-c(1,4,7,10)]
#DIse <- temp$SE[-c(1,4,7,10)]
temp<-read.table("compDImoments_agg.txt")
compDImoments_simple<-temp$EST[]
compDIse_simple<-temp$SE[]

temp<-read.table("compDImoments.txt")
compDImoments<-temp$EST[]
compDIse<-temp$SE[]
#compDImoments<-temp$EST[-c(5,6,11,12)]
#compDIse<-temp$SE[-c(5,6,11,12)]

setwd("../eventstudy/PSID_eventstudy")
rform_events_mod<-as.data.table(read.csv("./mod_eventmeans_table.csv"))
rform_events_sev<-as.data.table(read.csv("./sev_eventmeans_table.csv"))
eventDImoments_single<-c(unlist(rform_events_mod[X%in%c("as.factor(married)0"),'Estimate']),
                         unlist(rform_events_sev[X%in%c("as.factor(married)0"),'Estimate']))
eventDIse_single<-c(unlist(rform_events_mod[X%in%c("as.factor(married)0"),'Std..Error']),
                    unlist(rform_events_sev[X%in%c("as.factor(married)0"),'Std..Error']))
eventDImoments_married<-c(unlist(rform_events_mod[X%in%c("as.factor(married)1"),'Estimate']),
                         unlist(rform_events_sev[X%in%c("as.factor(married)1"),'Estimate']))
eventDIse_married<-c(unlist(rform_events_mod[X%in%c("as.factor(married)1"),'Std..Error']),
                    unlist(rform_events_sev[X%in%c("as.factor(married)1"),'Std..Error']))

eventDImoments<-c(eventDImoments_single,eventDImoments_married)
eventDIse<-c(eventDIse_single,eventDIse_married)

rform_events_mod<-as.data.table(read.csv("./mod_agg_eventmeans_table.csv"))
rform_events_sev<-as.data.table(read.csv("./sev_agg_eventmeans_table.csv"))
eventDImoments_simple<-c(unlist(rform_events_mod[X%in%c("Estimate"),'x']),
                         unlist(rform_events_sev[X%in%c("Estimate"),'x']))
eventDIse_simple<-c(unlist(rform_events_mod[X%in%c("Std. Error"),'x']),
                  unlist(rform_events_sev[X%in%c("Std. Error"),'x']))

setwd("../../output")

# temp1<-read.table("eventheadmodDImoments_ddbal.txt")
# temp2<-read.table("eventheadsevDImoments_ddbal.txt")
# temp1[,"V3"]<-as.character(temp1[,"V3"])
# temp1[,"V3"]<-as.numeric(gsub("\\(","",gsub("\\)","",temp1[,"V3"])))
# temp2[,"V3"]<-as.character(temp2[,"V3"])
# temp2[,"V3"]<-as.numeric(gsub("\\(","",gsub("\\)","",temp2[,"V3"])))
# eventDImoments<-c(temp1[1,'V2'],temp2[1,'V2'],temp1[2,'V2'],temp2[2,'V2'])
# eventDIse<-c(temp1[1,'V3'],temp2[1,'V3'],temp1[2,'V3'],temp2[2,'V3'])

temp<-read.table("flowDImoments_agg.txt")
flowDImoments_simple<-temp$EST[]
flowDIse_simple<- temp$SE[]


temp<-read.table("flowDImomentsbymar.txt")
flowDImoments<-temp$EST[]
flowDIse<- temp$SE[]
#flowDImoments<-temp$EST[-c(1,4)]
#flowDIse<- temp$SE[-c(1,4)]

temp<-read.table("contflowDImoments_agg.txt")
contflowDImoments_simple<-temp$EST
names(contflowDImoments_simple)<-rownames(temp)
contflowDIse_simple<- temp$SE
names(contflowDIse_simple)<-rownames(temp)

temp<-read.table("contflowDImomentsbymar.txt")
contflowDImoments<-temp$EST
names(contflowDImoments)<-rownames(temp)
contflowDIse<- temp$SE
names(contflowDIse)<-rownames(temp)


temp<-read.table("jointworkmoments.txt")
jointworkmoments<-temp$EST #(no work, no spwork), (no work, sp work), (work, no spwork), (work, spwork)

temp<-read.table("jointworkpartmoments.txt")
jointworkpartmoments<-temp$EST #(no work, no spwork), (no work, sp work), (work, no spwork), (work, spwork)
#Including event moments instead of flow and composition moments:
datamoments<-c(wageregmoments,spwageregmoments,wagevarmoments,
               consmoments_simple_all,workmoments_simple,spworkmoments, 
               #stockDImoments,
               compDImoments_simple,
               flowDImoments_simple,
#               contflowDImoments
               eventDImoments_simple
               )


#weightmat<-diag(c((1/(consse^2)),(1/(workse^2)), (1/(workpartse^2)),(1/(spworkse^2)), (1/(spworkpartse^2)),(1/(DIse^2)),(1/(compDIse^2)),(1/(flowDIse^2))))*1e8
##normalizing weightmat: consumption regression, work ,and DI moments each equally weighted
wagemat <- c(1/wageregse^2,1/spwageregse^2,1/wagevarse^2)
consmat <- (1/(consse_simple_all^2))#/sum((1/(consse_all^2)))
workmat <- c((1/(workse_simple^2)),(1/(spworkse^2)))
#workmat<-workmat /sum(workmat)
#Including event moments instead of flow and composition moments:
healthyflowse_simple<-1/(flowDIse_simple[grepl("L=0",names(flowDIse_simple))]^2)
healthyflowse_simple<-healthyflowse_simple/sum(healthyflowse_simple)
otherflowse_simple<-1/(flowDIse_simple[!grepl("L=0",names(flowDIse_simple))]^2)
otherflowse_simple<-otherflowse_simple/sum(otherflowse_simple)
DImat <- c(#(1/(DIse^2)),
           (1/(compDIse_simple^2)),
           (1/(flowDIse_simple^2)),
           #(1/(contflowDIse^2))
           (1/(eventDIse_simple^2))
            )
# DImat <-c(healthyflowse[1],
#           otherflowse[1:2],
#           healthyflowse[2],
#           otherflowse[3:4],
#           healthyflowse[3],
#           otherflowse[5:6],
#           healthyflowse[4],
#           otherflowse[7:8]
#           ) 
#DImat <- rep(1,length(DImat))
#DImat <- DImat/sum(DImat)
weightmat <- diag(c(wagemat,consmat,workmat,DImat))*1e8
numS=10
numsims<-read.table(paste0("numtypes",typename,".txt"))
numsims<-numsims[!rownames(numsims)=="Total",]*numS

weightmat[which(is.infinite(weightmat))]<-max(weightmat[which(!is.infinite(weightmat))])

#NOTE: need to adjust more than just this parameter if changing number of house types
numtypes<-length(numsims)
agelength<-1


#########################
#PARALLEL STUFF FOR LINUX:
#########################
particles<-mpi.universe.size()-1 # no. particles in the particle swarm search
print(mpi.universe.size())
#numThreads<-10
seed=6327
partic_mult<-1
##partic_mult<-4
##threads<-10
#########################




marriageprob_init<-as.data.table(read.csv(paste0("marriagemeans",typename,".csv"),row.names=NULL))
marriageprob_init<-marriageprob_init[,married]



states<-statesetup(dir=".",filesuffix=typename,
                   asset_gridsize=20,wage_gridsize=10)

#states<-statesetup(dir="/home/mdkellogg/StructuralModel",filesuffix=typename,
#                   asset_gridsize=20,wage_gridsize=10)
Asset<-states$Asset
foodstamps<-states$foodstamps
numkids<-states$numkids
healthprob<-states$healthprob
marriageprob<-states$marriageprob
for(aa in 1:3){
  for(bb in 1:3){
    marriageprob[[aa]][[bb]][,1]<-1
    marriageprob[[aa]][[bb]][,2]<-0
  }
}#marriageprob_init<-states$marriageprob_init
oecd<-states$oecd

#setwd("../../StructuralModel")
setwd("../modeloutput")

# Asset<-array(NA,c(numtypes,50*agelength+1,asset_gridsize))
# for (typ in 1:numtypes){
#   assetlow<-rep(0,50*agelength+1)
#   assethigh=NULL
#   assethigh[1] <- 2000
#   Asset[typ,1,]<-seq(from=assetlow[1],to=assethigh[1],length.out=asset_gridsize)
#   for(t in 2:40*agelength){
#     assethigh[t] <- max(assethigh[t-1]+2000*max(Wage[typ,t,])/(t^0.5),2000)
#     Asset[typ,t,]<-seq(from=assetlow[t],to=assethigh[t],length.out=asset_gridsize)
#   }
#   for(t in (40*agelength+1):50*agelength){
#     assethigh[t] <- max(assethigh[t-1]-2000*max(Wage[typ,t,])/(t^0.5),2000)
#     Asset[typ,t,]<-seq(from=assetlow[t],to=assethigh[t],length.out=asset_gridsize)
#   }
#   
#   Asset[typ,50*agelength+1,]<-Asset[typ,50*agelength,]
# }

# Ctime<-NULL
# Ctime[1] <-Sys.time()
# output<-Welfare_solve(as.vector(initguess),datamoments[,1],invvarcov,healthprob,foodstamps,numkids,Asset,0,1)
# Ctime[2] <-Sys.time()
# Rtime<-NULL
# Rtime[1] <-Sys.time()
# outputR<-welfare(initguess,datamoments,invvarcov,healthprob,foodstamps,numkids,Asset)
# Rtime[2] <-Sys.time()

#########################
#PARALLEL STUFF FOR LINUX:
#########################
cl <- makeCluster(particles, type="MPI",outfile="")
registerDoSNOW(cl)
clusterEvalQ(cl,{
  library(Rcpp)
  library(data.table)
  library(RcppArmadillo)
  library(RcppParallel)
  library(ggplot2)
  library(parallel)
  library(snow)
  library(doSNOW)
  library(Rmpi)
  library(spousalVFIParallel)
  library(unixtools)
  library(readr)
  library(fixest)
  library(readstata13)

  source("Welfare_ParallelSpouseC.R")}
)

#########################

#(rough) result with Voena-type preferences:
#"Objective evaluated at 1.54310315225894"
roughresult<-c(theta=-0.422571932993375,sptheta=-0.477994837409422,eta=-0.333022576732923,speta=0.111072538682675,
               etaspeta=0,delta=0.254196342081335,F1=0.625234518499814,F2=0.941136111049607,spF2=1.07713807603117,pi0.young=0.00599999709566033,
               pi1.young=0.170970905428669,pi2.young=0.320558946872059,pi0.old=0.0754553925296934,pi1.old=0.205705982189693,pi2.old=0.454573696903444)

initialguess<-c(theta=-0.448, sptheta=-0.448,eta=c(-0.185), speta=c(-0.185),
                 delta=0.062
                ,F1=0,F2=0.547,F3=0.952,spF1=0.25, spF2=0.952,
                pi0.young=0.006,pi1.young=0.171,pi2.young=0.331,pi0.old=0.075,
                pi1.old=0.180,pi2.old=0.626)

# newguess<- c(-0.467776561,	-0.030043986,	-0.003019146,	-0.625170263,
#              -0.100576715,	1.370785275,	0.063538949,	0.002756961,
#              0.502114648,	2.749056557,	1.579133453,	1.952137472,
#              0.625790516,	0.518175796,	0.804022976,	0.647171689,	0.599334807,	0.984376506)

#Guess, killing flexibility in theta
# newguess<- c(-1,	-0.1,	-0.086249373,	-0.086249373/2,	-0.075519018, -0.075519018/2,
#   0.087225554,	0,	0.446291087,	1.143357757, 0, 0.446291087/2, 1.143357757/2,
#   1.083249787, 1.223382484, 1.083249787/2, 1.223382484/2,	0.286825794,	0.338727586,	0.613171924,
#   0.394304924,	0.409813551,	0.83582787)
typename1<-"stationary"
fullpar<-NULL
for(particle in 1:(5)){
  guessvals<-read.csv(paste0('./particlefiles/outputvals_fixmar',typename1,particle,'.csv'),header=FALSE)
  guessvals<-apply(guessvals,MARGIN=2, function(x) as.numeric(as.character(x)))
  guessvals<-guessvals[,-1]
  guessvals<-guessvals[guessvals[,ncol(guessvals)]==999,]
  guessvals<-guessvals[,-ncol(guessvals)]
  newguess<-unlist(as.vector(guessvals[which.min(guessvals[,dim(guessvals)[2]]),]))
  fullpar<-cbind(fullpar,newguess)
}
# guessvals<-read.csv(paste0('./particlefiles/outputvals_short2_',typename1,'.csv'),header=FALSE)
# guessvals<-apply(guessvals,MARGIN=2, function(x) as.numeric(as.character(x)))
# guessvals<-guessvals[,-1]
# guessvals<-guessvals[guessvals[,ncol(guessvals)]==999,]
# guessvals<-guessvals[,-ncol(guessvals)]
newguess<-unlist(as.vector(guessvals[which.min(guessvals[,dim(guessvals)[2]]),]))
fullpar<-cbind(fullpar,newguess)

newguess<-fullpar[-dim(fullpar)[1],which.min(fullpar[dim(fullpar)[1],])]
fullpar<-fullpar[-dim(fullpar)[1],]
#fullpar[,2]<-newguess
#fullpar[(nrow(fullpar)-11):nrow(fullpar),2]<-0
#adding in work probabilities:
#newguess<-c(newguess,newguess[(length(newguess)-11):length(newguess)])

newguess<-c(#0,-0.15,0.074,-0.074,0.895,0.889,1.516,
            #0.02,0.103,-0.084, 0.550,0.511,0.824,
            #0.142^2,0.1807^2,0.203,0.202^2,0.227^2,
            newguess)
#newguess<-c(newguess[1:7],0,0,newguess[8:9],0,newguess[10:length(newguess)])
names(newguess)<-c("mod","sev","age","age2","Type1","Type2","Type3",
                   "spmod","spage","spage2","spType1","spType2","spType3",
                   "wagevar","spwagevar","wagecorr","noisevar","spnoisevar",
                   "theta_mod","theta_sev",
                   "spmu","sptheta",
                   "eta","etapart","eta_single","etapart_single","eta_mod","eta_sev",
                   "speta","spetapart","speta_mod",
                   "delta","delta_single", "spdelta",
                   "F1","F2","F3","Fpart1","Fpart2","Fpart3","Fold","Fpartold","F1_single","F2_single","F3_single","Fpart1_single","Fpart2_single",
                   "Fpart3_single","Fold_single","Fpartold_single","spF1","spF2","spFpart1","spFpart2","spFold","spFpartold",
                   "pi0.young.married","pi1.young.married","pi2.young.married","pi0.old.married","pi1.old.married","pi2.old.married",
                   "pi0.young.single","pi1.young.single","pi2.young.single","pi0.old.single","pi1.old.single","pi2.old.single",
                   "pi0.young.married.work","pi1.young.married.work","pi2.young.married.work","pi0.old.married.work","pi1.old.married.work","pi2.old.married.work",
                   "pi0.young.single.work","pi1.young.single.work","pi2.young.single.work","pi0.old.single.work","pi1.old.single.work","pi2.old.single.work")

#rownames(fullpar)<-names(newguess)
#fullpar["Fpartold",]<-seq(from=0,to=0.5,length.out=partic_mult*particles)
#newguess<-fullpar[,1]

upper<-c(0.2,0.2,0.2,0.2,3,3,5, #Ref person work params
         0.05,0.2,0.2,1,1,2, #spouse work params
         c(1,1,0.5,1,1), #work variance params
         0,0, #Thetas: disutility of health
         1.5,0.2, #spmu and sptheta: disutility of marriage and spousal health
         0,0, #etas: disutility of work, married
         0,0, #etas: disutility of work, single
         0, 0, #etas: disutility of work, when disabled
         0.2,0,0, #etas: disutility of work, for spouse
         0.2,0.6,1, #deltas: job loss probabilities
         0,1,3, #F: fixed cost of full-time work, married
         0,0,0, #F: fixed cost of part-time work, married
         #         0.5,2,4, #F: fixed cost of part-time work, married
         0.5,0, #F: extra fixed cost of work for being old, married
         #0.5,2, #F: extra fixed cost of work for being old, married
         0,2,3, #F: fixed cost of full-time work, single
         0,0,0, #F: fixed cost of part-time work, single
         #         0.5,2,3, #F: fixed cost of part-time work, single
         0.5,0, #F: extra fixed cost of work for being old, single
#         0.5,2, #F: extra fixed cost of work for being old, single
         0,2, #F: fixed costs of full-timework for spouse
         0,0, #F: fixed costs of part-timework for spouse
         #         3,3, #F: fixed costs of part-timework for spouse
         0.5,0, #F: extra fxed cost of work for being old, spouse
        #0.5,1, #F: extra fxed cost of work for being old, spouse
        0.4,1,1,1,1,1, #allowance probs, married
        0,0,0,0,0,0, #allowance probs, single
#         0.1,0.1,0.1,0.1,0.1,0.1, #allowance probs, married working
#         0.1,0.1,0.1,0.1,0.1,0.1) #allowance probs, single working
         0,0,0,0,0,0, #allowance probs, married working
         0,0,0,0,0,0 #allowance probs, single working
) 

lower<-c(rep(-1,7), #Ref person work params
         rep(-1,6), #spouse work params
         rep(0.0001,5), #work variance params
         -1.5,-1.5,
        0.2,-1,
#         -0.5,-0.5,
         -0.5,0,
#            -1,-0.5,
            -1,0,
            -0.2,-1,
          #-1,-1,-1,
          -1,0,-1,
         0,0,0,
         0,0,0,
         0,0,0,
         0,0,
         0,0,0,
         0,0,0,
         0,0,
         0,0,
         0,0,
         0,0,
         0,0,0.1,0,0,0.3,
         0,0,0,0,0,0.1,
         0,0,0,0,0,0,
         0,0,0,0,0,0
)


#dropping spousal parameters:
#upper[c(3,10:12,15,32:37)]<-newguess[c(3,10:12,15,32:37)]
#lower[c(3,10:12,15,32:37)]<-newguess[c(3,10:12,15,32:37)]

#upper<-newguess
#lower<-newguess
#upper[(length(newguess)-11):length(newguess)]<-1
#lower[(length(newguess)-11):length(newguess)]<-0
#upper[c(1,2,3,4,6,8,9,10,12)] <- 0
#lower[c(1,2,3,4,6,8,9,10,12)] <- -1
names(upper)<-names(newguess)
names(lower)<-names(newguess)

#Shutting downs spousal PT parameters
newguess[c("spetapart","spFpart1","spFpart2","spFpartold")]<-0
upper[c("spetapart","spFpart1","spFpart2","spFpartold")]<-0
lower[c("spetapart","spFpart1","spFpart2","spFpartold")]<-0

#Shutting down moderate theta and eta:
newguess[c("theta_mod","eta_mod", "spmu", "eta_single", "delta_single", "F2", "Fold", "F2_single", "F3_single", "Fold_single", "spFold",
           "pi0.young.single","pi1.young.single","pi2.young.single","pi0.old.single","pi1.old.single","pi2.old.single")]<-0
lower[c("theta_mod","eta_mod", "spmu", "eta_single", "delta_single", "F2", "Fold", "F2_single", "F3_single", "Fold_single", "spFold",
        "pi0.young.single","pi1.young.single","pi2.young.single","pi0.old.single","pi1.old.single","pi2.old.single")]<-0
upper[c("theta_mod","eta_mod", "spmu", "eta_single", "delta_single", "F2", "Fold", "F2_single", "F3_single", "Fold_single", "spFold",
        "pi0.young.single","pi1.young.single","pi2.young.single","pi0.old.single","pi1.old.single","pi2.old.single")]<-0



##shutting down allowance prob for singles, setting it same as married
#newguess[(length(newguess)-5):length(newguess)]<-0
#lower[(length(newguess)-5):length(newguess)]<-0
#upper[(length(newguess)-5):length(newguess)]<-0


newguess[which(newguess>upper)]<-upper[which(newguess>upper)]
newguess[which(newguess<lower)]<-lower[which(newguess<lower)]



# #shutting down everything but allowance probs:
# lower[1:37]<-newguess[1:37]
# upper[1:37]<-newguess[1:37]

##Shutting down PT work moments:
#newguess[grepl("part",names(newguess))]<-0
#upper[grepl("part",names(newguess))]<-0
#lower[grepl("part",names(newguess))]<-0


# bwidth<-0.05
# upper<-pmin(upper,newguess+bwidth)
# lower<-pmax(lower,newguess-bwidth)

#fullpar<-sapply(1:particles,FUN=function(x) newguess)
#rownames(fullpar)<-names(newguess)
#fullpar["sptheta",]<-seq(from=newguess["sptheta"],to=lower["sptheta"],length.out=particles)
#fullpar["spdelta",]<-seq(from=newguess["spdelta"],to=upper["spdelta"],length.out=particles)


set.seed(seed)
.Random.seed<-c(10407L,  -888780854L,   762756803L,   780773672L,  -352617207L,  1841665334L, -2063451585L)
#partnershocks<-exp(rmvnorm(50*sum(numsims),c(0,0),sigma=wagecov))
partnershocks<-rmvnorm(50*sum(numsims),c(0,0))
shocks.wageraw <- partnershocks[,1]
shocks.spwageraw <- partnershocks[,2]
shocks.wagenoise <- rnorm(50*sum(numsims),0,1)
shocks.spwagenoise <- rnorm(50*sum(numsims),0,1)
shocks.health<- runif(50*sum(numsims))
shocks.allow<- runif(50*sum(numsims))
shocks.reassess<- runif(50*sum(numsims))
shocks.jobloss<- runif(50*sum(numsims))
shocks.spjobloss<- runif(50*sum(numsims))
shocks.marriage<- runif(50*sum(numsims))


select<-drawselect(numsims)
selectC<-drawselectC(numsims)

 setThreadOptions(numThreads=10,stackSize=1e07)

# CLEAR THE OUTPUT FILE IN ADVANCE---NOT FOR NOW.
#if(file.exists("outputvals.csv")) file.remove("outputvals.csv")
#########################
#PARALLEL STUFF FOR LINUX:
#########################

ex<- Filter(function(x) is.function(get(x,.GlobalEnv)),ls(.GlobalEnv))
clusterExport(cl,ex)

clusterEvalQ(cl,{
  setThreadOptions(numThreads=10,stackSize=1e07)
}
)

#########################

result<- foreach(i = 1:(partic_mult*(particles+1)))%dopar%{
  psoptim_subparallel(newguess,welfareC,
#              psoptim(newguess,welfareC,
##                      fullpar = fullpar,
                      lower=lower,
                          upper=upper,
                      outfile=paste0('./particlefiles/outputvals_fixmar',typename),
                      particle=i,
                      specification='base', parnames=names(newguess), datamoments=datamoments,
                          healthprob=healthprob,marriageprob=marriageprob,foodstamps=foodstamps,numkids=numkids,
                          weightmat=weightmat,Asset=Asset,Wage=Wage, spWage=spWage,agelength=agelength,
                           wagenoise=shocks.wagenoise, spwagenoise=shocks.spwagenoise,
                          shocks.wage=shocks.wageraw,shocks.spwage=shocks.spwageraw, shocks.health=shocks.health, shocks.allow=shocks.allow,
                         shocks.reassess=shocks.reassess, shocks.jobloss=shocks.jobloss, shocks.spjobloss=shocks.spjobloss,
                         shocks.marriage=shocks.marriage, marriageprob_init=marriageprob_init, numsims=numsims, lumpsum = 0,
                         select=select,selectC=selectC, seed=seed, 
                         control=list(s=partic_mult*particles,maxit=1000),
                         noflow=FALSE,nocomp=FALSE,nostockDI=TRUE,noeventDI=FALSE,nocontflow=TRUE,
                         workDI=0,noSSDI=0,onlyFT=1,
                         single_sameallowanceprob=TRUE,
                         work_noallowanceprob=FALSE,
                         single_sameprefs = TRUE,
                         health_simpleprefs = TRUE,
                         nospmoments=FALSE,
                         noworkDImoments=TRUE,
                         consmoments_simple=TRUE,workmoments_simple=TRUE,DImoments_simple=TRUE,
                         onlymarried=FALSE, onlysingle=FALSE)
}



stopCluster(cl)

